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ABSTRACT 

As is common for antenna arrays in radio astronomy, the output of the MWA's correlator is the 
intensity measured in visibility space. In addition, the final power spectrum will be created in visi- 
bility space. As such, correcting for the ionosphere in visibility space instead of real space saves the 
computation required to inverse Fourier transform to real space and then Fourier transform back (a 
significant decrease in computation for systems operating in real time such as the MWA.) In this paper, 
we explore this problem of correcting for ionospheric distortions in the uv-plane. The mathematical 
formula for obtaining the unperturbed data from that reflected by the ionosphere is non-local, which 
in any practical application creates edge effects because of the finite nature of the wi>-plane (section 
1.2). In addition, obtaining an analytic solution for the unperturbed intensity is quite difficult, and 
can only be done using very specific expansions of ionospheric perturbations. We choose one of these 
models (with perturbations as sinusoidal modes, section 2) and run numerical codes to further study 
the correction. Numerically implementing this correction to too few orders distorts the data in such 
a way as to be worse than not correcting at all (section 4). It is therefore critical to correct to the 
correct number of orders, and we present an analytic estimate for the optimal order (section 5). This 
analytic estimate shows that the optimal number of orders varies with u, and in particular increases 
as u increases along the direction of an ionospheric distorting mode. Based on this observation, we 
then investigate a couple of methods which save computation (section 6). These methods are (a) 
eliminating the intensity at values of u which require too many orders, and (b) correcting to different 
orders at different u. Both methods prove successful, although the first creates a loss of some precision 
in the real space sky. We conclude by considering an alternate form with which to model ionospheric 
perturbations (section 7). This alternate form was once again chosen because it lends itself to an 
analytic solution, but contains as many (if not more) downfalls than the original choice. 
Subject headings: 



1. BACKGROUND ON THE IONOSPHERIC CORRECTION 
OPERATOR 

We begin with the basic background on the iono- 
spheric operator, followed by that of the ionospheric 
correction operator. These sections expand upon the 
mathematical framework created and briefly outlined in 
Morales, et. al., (1). 

1.1. Ionospheric Operator 

The ionospheric operator A(8',9) is the operator 

which takes an unperturbed map of the sky 1(9) and 

maps it to a perturbed map 1(9') which has been dis- 
torted by the ionosphere, 



here, and with other operators later, is indicative of the 
direction of change.) 

In the regime of MWA, it is a very good approximation 
that the mapping of angles is approximately linear with 
only a small deviation 69(9, t), 



69(8, t). 



(1.2) 



1(9') = A(9', 9)1(9). 



1.1 



Put simply, A(9',9) is a generalized coordinate change 
from 6 to 9'. (The order of the arguments of A(9',9) 
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The perturbed intensity 1(9') at 9' is the summed con- 
tribution from the intensities I(9i) at all 9i where this 
relation holds; ie, 

1$) = E ^iW' - £ - (1-3) 



where 5 represents the Dirac delta function. For example, 
if (only) 8\ and 92 are mapped to 9[ according to equa- 
tion 1.2, then I(9\) = I(9[) + J(#l). In the limit that 
9 is a continuous variable, this discrete sum for 1(9') is 



2 



altered to an integral, 

1(6') = J d 2 0A(0', 0)1(0), (1.4) 

where 

A(6',8) = 8(8' - 8- 80(6, t)). (1.5) 

(Notice that we've used that the magnitude of the delta- 
function's argument's gradient is approximately 1 here). 
From this equation it is evident that in the limit of con- 
tinuous 9, the ionospheric operator equation becomes 

1(0') = A(0' ,0)1(0) = f d 2 0A(0', 0)1(0). (1.6) 



So far our calculations have been confined to real space, 
using the variables { 0' , }. However, as is common for 
antenna arrays in radio astronomy, the output of the cor- 
relator for MWA will actually be the Fourier transform 
of the real space sky intensity, I(u). In addition, the final 
power spectrum of the sky will also be measured in this 
Fourier transfrom space (also known as visibility space, 
or uv -space). Therefore, correcting for the ionosphere in 
real space requires inverse Fourier transforming to real 
space, making the correction, and then Fourier trans- 
forming back to visibility space. The problem with this 
is that Fourier transforming is computationally expen- 
sive, especially for a system operating in real time, such 
as the MWA's Real Time System (4). As such, correcting 
for the ionosphere in the uv-plane would greatly reduce 
computation. We now study the nature of this uw-plane 
correction. 

Let I(u) represent the Fourier transform of 1(6). De- 
fine this Fourier transform by 

1(5) = F(u,0)I(0) = J d 2 0e- tZ - S I(0), (1.7) 

and the corresponding inverse Fourier transform by 

d 2 u 



1(0) = F~ (6,u)I(u) 



(2nJ- 



;e lv -"I(u). (1.8) 



Now define A(u',u) to be the ionospheric operator in 
the uu— plane; that is, the operator that maps the un- 
perturbed map I(u) to the perturbed map I(u'), 



I(u') = A(u',u)I(u). 



(1.9) 



One way to obtain I(u') from I(u) is to inverse Fourier 
transform I(u) to 1(6) using F^ 1 (8,u), apply the iono- 
spheric operator A(6',0) in real space to obtain 1(6'), 
and then Fourier transform 1(6') to I(u') using F(u', 6'). 
In all, 

J(u') ^F(u',0')A(0',0)F- 1 (0,u)I(u). (1.10) 

Comparison of this to the definition of A(u',u) shows 
that 

A(u',u) = F(u',0')A(0',0)F- 1 (0,u). (1.11) 

This expression may be thought of as simply a basis 
change of A from {0,0'} to {u, u'}. The three opera- 
tors on the right here have all been previously given. 



Plugging in these predetermined expressions (equations 
1.6, 1.7, and 1.8) and simplifying as much as possible, we 
obtain 



(1.12) 

d 2 9e i(u-u')-e-iu'-se \ /(£)(i. i3) 



I(u')=A(u',u)I(u) 

= /(2^(/ 

Notice that the integral over 0' has been evaluated by 
using the delta function from the expression for A(0', 0) 
(see equations 1.5 and 1.6; here again we use that the 
magnitude of the delta-function's argument's gradient is 
approximately 1). One interesting characteristic of this 
expression is the non-local nature, by which finding the 

value of the perturbed intensity I(u') at one particular 
u' requires knowing the value of the pure intensity I(u) 
at other u ^ u' . This property will also appear in the 
ionospheric correction operator, found below. 

1.2. Ionospheric Correction Operator 

The previous section corresponds to the operator 
which distorts the pure data into the perturbed data, but 
the reverse process is what actually interests us - we want 
to correct for the effect of the ionosphere to obtain the 
pure data from the perturbed data. Define A T (0,6') to 
be this ionospheric correction operator which corrects for 
the influence of the ionosphere by mapping the perturbed 
map of the sky 1(8') back to the unperturbed map 1(6), 

I(0) = A T (0,0')I(0'). (1.14) 

The MWA will not run during periods of scintillation (at 
which times multiple values of are perturbed to the 
same 0'), but will instead run during times when it is a 
very good approximation that the mapping from 0' to 
is one-to-one and approximately linear with only a small 
correction, 

8 = 8' + 88' (8', t). (1.15) 

The derivation of an expression for the ionospheric cor- 
rection operator in the itw-plane A T (u, u') follows analo- 
gously to the derivation of A(u', u), so I'll merely quote 
the result: 

I(u)=A T (u,u')I(u') 



i{u' —u)-0' — iu-< 



(1.16) 
/(«')■ (1-17) 



This is the most general expression for obtaining I(u) 

from I(u'); proceeding further requires knowledge of the 

ionospheric perturbation 88' . Numerically solving for 
I(u) using this equation is an incredibly daunting task 

for an arbitrary choice of 86' , as it involves a double inte- 
gral over all space for every value of u. Therefore, unless 

we find a choice of 86' which offers an analytic solution 
for I(u), the correction for the ionosphere in the uw-plane 
will actually be more computationally expensive than the 
two Fourier transforms necessary to correct for the iono- 
sphere in real space. Unfortunately, choices of 86' which 
lend themselves to analytic solutions are hard to come by. 
There are a couple, however, and they will be discussed 
in the following sections. 
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2. A SPECIFIC FORM FOR SB': SUM OVER SINUSOIDAL 
MODES 

The above equation for I{u) contains an exponen- 
tial with 59' in the exponent. By expanding this expo- 
nential, we obtain a form for I(u) which may be solved 
analytically for a couple of choices of 59'. To be explicit, 
expanding the exponential in 59' , 



' = E 



(-iu-59') r 



n=0 



ni 



(2.1) 



leads to 



(Interchanging an infinite sum and an integral requires 
that the sum be uniformly convergent, which will be true 
for all choices of 59' that we choose.) 

With the above expansion, I(u) may be solved analyt- 
ically if we choose 



M 



5fr = Re^iaJLe-*™-*'). 



(2.3) 



m— 1 



Here, the b m are chosen to be purely real, but the (Xyyi are 
allowed to assume complex values. Physically, this choice 

of 59' corresponds to modeling the integral along the line 
of sight of the ionosphere's electron density n e (9',h) (h 
is the distance along the line of sight) as a sum over 
sinusoidal modes, 



/ 



M 



(9',h)dh oc J2 



-ib m -9< 



(2.4) 



The reflection by the ionosphere is then related to this 
by 



59' oc V ffl (^j n e (9',h)dh 



M 

E 

771—1 



-ibr, 



-ib„ 



(2.5) 



where V e -, represents the two-dimensional gradient with 
respect to (9' x , 0' ). (The actual shift, of course, is a Her- 
mitian observable, so only the real part of this is included 

in 59 .) Notice that this choice ultimately stems from 
our decision to model density fluctuations from the iono- 
sphere as sinusoidal modes. Actually, any orthonormal 
basis would have sufficed here. Once again, this particu- 
lar choice was made because it allows an analytic solution 
for the unperturbed intensity I{u). (An alternate choice 
which likewise offers an analytic solution will be briefly 
discussed later on in section 7.) 

With this choice of 89' , the intensity I(u) becomes 



v n- 

I(u'). (2.7) 



■Re(^2 iambme- ih ™- e ') 

ra—l 



The math leading to a solution for I(u) may be found 
in the appendix, section A. Very briefly, the integral is 
solved by conveniently redefining the ionospheric modes 
{a m , b m } (as given below), performing a multinomial ex- 
pansion on the term raised to the power n, recognizing 

that the final product of this expansion leaves the 9' in- 
tegral in the form a delta function, and then using that 
delta function to evaluate the integral over u' . The final 
solution is 

oc / (1M , _ r w \ 

'(*) = E E n^ ^+EU). 



n=0 h,h, ---thM \<?=1 



where 



c<,= ? 



\a q if q < M + 1 



ia * q _ M iiq>M + l 



d q = 



?q-M 



ifq<M + l 
\{q>M+l, 



(2.8) 
(2.9) 

(2.10) 



and the summation over all l q (with 2M > q > 1) is a 
restricted sum such that J2q K — n an( i ^ > 0. This 
equation might look a little daunting, but it may be 
thought of simply as the addition of many delta func- 
tions of varying amplitudes, with those further from the 
point u in question tending to contribute less to the sum. 
(This is, in fact, similar to what one sees with intermod- 
ulation distortion). Notice that although there are M 
modes distorting the sky, the sums and products above 
involve 2M > q > 1. Thus, there appear to be 2M ef- 
fective modes distorting the sky. This factor of 2 comes 
from the constraint that 89' be real, as may be more eas- 
ily seen by following the math provided in the appendix, 
section A. 

Another important feature of this solution is that it 
is inherently non-local, with the corrected intensity I(u) 
at a given u depending on the values of the perturbed 
intensity I(u') at the appropriate neighboring points u' — 
u + Iqdq- This non-locality, which is also evident in 
the most general form for I(u) (equation 1.17), will create 
edge effects because the uw-plane is finite in all practical 
applications, as will be more easily seen and understood 
later in section 4.2. 

It should also be pointed out that this equation for I(u) 
is the result of a double expansion: the Taylor-series ex- 
pansion of the exponential containing 89' (see equation 
2.1), which is now evident in the summation over n, and 

the expansion of the ionospheric perturbation 59' itself 
into sinusoidal modes (see equation 2.3), which is now 
evident in the restricted sum over l q . Throughout this 
paper we will assume that this second expansion is "per- 
fect" ; that is, we will assume that we are able to perfectly 
model the ionosphere with the M modes that we assume 
are provided for us. We will instead study the errors 
created by truncating the first expansion. 

3. POTENTIAL PROBLEM: COMPUTATIONAL 
FEASIBILITY OF THE TWO EXPANSIONS 

As previously stated, the main goal of the uu-planc 
correction is to correct for the ionosphere in a less com- 
putationally intensive manner than that required for the 
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real space correction. Our final expression for I(u), how- 
ever, contains an infinite sum over n. Clearly, making the 
uu-plane correction computationally feasible will require 
truncating this sum after a finite number of terms. The 
next section will explore the effect of such a truncation. 

But even truncating this sum over n cannot guarantee 
the computational feasibility of the icu-plane correction 
because of the second expansion over sinusodial modes 
and its resulting restricted summation. More specifically, 
the number of terms in the restricted sum over all pos- 
sible combinations of l q such that Y] q l q — n may be 
calculated through the following trick: if 2M represents 
the total number of effective modes, then consider the 
problem of arranging n balls and 2M — 1 partitions in a 
straight line. Here the i-th partition marks the stopping 
point where U ends and begins. For example, if for 
a particular arrangement 7 balls lie between the 4th and 
5th partitions, then l 5 — 7 for that arrangement. The to- 
tal number of ways to arrange these n + 2M — 1 objects is 
(n + 2M — 1)1. Of course, exchanging the position of any 
two of the same object (ball or partition) does not lead 
to a different arrangement, so the total number of terms 



in the restricted sum such that ^ 



2AI 

# of terms in the sum l q 

q = l 



2M 



1 



n is given by 

(n + 2M- 1)! 
n\{2M-l)[ ■ 

(3-1) 



As an example, suppose that we wish to calculate this 
sum for 10 modes (2M = 20 effective modes) to the 40th 
order in n. Using the above formula, we calculate that 
such a sum has approximately 70 trillion terms. From 
this we see that the uu-plane correction is only compu- 
tationally feasible if the number of modes necessary to 
model the sky M and the number of orders necessary in 
the expansion of the exponential n are relatively small. 

4. COMPUTATIONAL RESULTS: TRUNCATING THE 
INFINITE SUM OVER N 

As stated in the previous section, the wu-plane cor- 
rection is only computationally feasible if we truncate 
the infinite sum over n. At this point we pause to study 
the results of truncating this sum after a finite order of 
correction, n = n max . 

4.1. Simple Sky Model 

In order to study the qualitative effects of truncat- 
ing the sum over n, we created a simple sky and per- 
turbed it with a simple mode, and then used numerical 
code to correct for this (known) perturbation in the uv- 
plane by using the mathematical formula found above in 
equation 2.8. 

The initial, unperturbed sky is shown in the top panel 
of figure 1. We will refer to this sky throughout this 
paper as the simple sky. It contains a 101 x 101 pixel 
array with a spacing of 4 arcmin between pixels, which is 
the approximate resolution we expect for the final MWA 
array (The axes in this and all the following real space 
figures are labeled in radians). This pure sky is a single 
source sky: The value of the intensity at all pixels is 
set to 0.0 except at one pixel where the value 1.0. (The 
important qualitative results found below would not be 
altered by including side lobes, so we will leave them out 
to keep things simpler.) 
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Fig. 1. — Top Panel A simple model of a sky (referred to as 
simple sky throughout this paper) consisting of a single source; ie, 
the value at each pixel is set to 0.0 except at one pixel, where it is 
set to 1.0. Notice that this simple model has no side lobes. Bottom 
Panel The simple sky from the top panel, but shifted by a single 
ionospheric mode, as given by equation 2.3 with a m = 1.0 X 10 -5 
radians squared and b m with a magnitude of 378.0 inverse radians 
and oriented in the 7r/4 direction with respect to the 9 X axis. It is 
difficult to see from these plots, but the single source has shifted in 
the 7r/4 direction by a few pixels and still has an intensity of 1.0. 
The axes of both graphs are measured in radians. 



Figure 2 shows the real part of the intensity in the 
uu-plane for this pure sky plotted in the third direc- 
tion (which is determined by the color scale) and demon- 
strates that the uncorrupted, wu-sky is a simple sinusoid 
(as one would expect for the Fourier-transform of a delta 
function). In contrast to this plot, the color scale for all 
the following uu-spacc plots is representative of the mag- 
nitude of the intensity I(u) in the uv-plane (although 
I(u) is complex, the important uu-plane results found 
below do not require phase information to understand 
graphically). The important feature to take away from 
this plot is that the absolute magnitude of the intensity 
is constant and of the order 10 -6 at all points in the uv- 
plane (although we've only plotted the real part here). 
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Fig. 2. — This is the pure visibility sky; ie, the Fourier transform 
of the simple unperturbed sky shown in the top panel in figure 
f. The dimensions of the axes are inverse radians, and the third 
dimension represented by the color scale is the real part of the 
intensity. Notice from these two panels that the pure, uncorrupted 
intensity in visibility space is a simple sinusoid, which is what one 
would expect for the Fourier transform of a delta-function. The 
important point to take away from this plot is that the magnitude 
intensity is on the order of 10 — 6 for our simple sky, and is the same 
for each uv-pixel (although only the real part is plotted here). All 
other ti?;-plots in this paper will have the magnitude of the intensity 
(and not just the real part) for the third, colored dimension. 

We then perturb this simple sky with a rather strong 
mode, as shown in the bottom panel of figure 1. This dis- 
torting mode has \a m \ — 1.0 x 1CT 5 radians squared and 

b m with a magnitude of 378.0 inverse radians and ori- 
ented in the 7r/4 direction as measured from the 9 X axis. 
(These values were chosen for the sole reason that they 
produce a strong shift of a few pixels (tens of arcmin- 
utes) , and thus accentuate the qualitative features of the 
ionospheric correction as seen below. A more realistic 
distortion will be discussed later in section 8.) Notice 
that the intensity is still 1.0 at exactly one pixel and 0.0 
at every other pixel, but now the location of this pixel 
has slightly shifted in the 7r/4 direction. (It should be 
noted that this ionospheric shift was applied to the con- 
tinuous sky with a delta function at one point, and not 
to its pixelized represention shown in the top panel of 
figure 1). Although we have not included the plot, the 
magnitude of the intensity in the perturbed ww-plane re- 
mains of the order of 10 -6 , its value in the pure visibility 
sky. 

4.2. Truncation Through 10 Orders for our Simple Sky 

We now begin correcting for this simple one-mode 
distortion using various values of n max . We (quite nat- 
urally) begin with the first order correction, n max = 1 
(n ma x = leads to no correction, see equation 2.8). Af- 
ter correcting to first order in the uw-plane, we inverse 
Fourier transform back to the real space sky shown in 
the top panel of figure 3 in order to determine the ef- 
fect that this first order correction has had on the real 
sky (in particular, we would like to know whether it has 
successfully shifted the single source back to its unper- 
turbed location). As it turns out, the correction to one 
order has not shifted the star from its perturbed location. 



The cross-like pattern of the star is somewhat interest- 
ing, but what is most important about this figure is that 
the maximum in the intensity has now doubled from 1.0 
to 2.0. 

A potential clue to this behavior is found by studying 
the visibility space sky corrected to first order, as shown 
in the bottom panel of figure 3. From this figure we see 
that the first order correction in the ww-plane has created 
an increase of an order of magnitude in the absolute value 
of the intensity at those points in the ww-plane furthest 
from the origin. 

Another important feature of this figure (although, as 
it turns out, it is not the cause of the increase in the real 
space intensity) is the ring around the outside edge of the 
figure. This ring is caused by the previously mentioned 
fact that the correction in the Mf-plane is non-local (see 
equation 2.8), combined with the finite nature of our nu- 
merical uu-plane. More specifically, points near the edge 
of our uu-plane may not obtain the full correction in vis- 
ibility space, because doing so requires pulling values of 
the intensity that are off the edge of the grid. Therefore, 
values near the edge are never fully corrected. We will 
later see that the result of this is a small spreading of the 
initial source (ie, a loss of precision), in real space. 

We now continue on to the second order correction 
(fi mol = 2). After inverse Fourier transforming, we ob- 
tain the real space sky shown in figure 4 (top panel). 
From this figure we see that the maximum in the inten- 
sity has increased even more, and is now 12 times its 
unperturbed and uncorrected value. The visibility sky 
after two orders of correction, as shown in the bottom 
panel of figure 4, has a maximum in the intensity that is 
now 100 times the value of the pure uv-sky. 

Correcting to higher orders, we see that the problem 
with intensities that are too high not only persists, but 
continues to get worse. After 5 orders, the intensity in the 
real space sky (figure 5, top panel) is 120 times too large, 
and that in the visibility space sky (figure 5, bottom 
panel) is 10 3 times too large; after 10 orders, the intensity 
in the real space sky (figure 6, top panel) is 600 times 
too large, and that in the visibility space sky (figure 6, 
bottom panel) is 10 4 times too large. In addition to 
the increase in the maximum in real space intensity, the 
source is also beginning to spread out and look less like 
a single point source. 

4.3. Making Sense of the Bizarre Behavior of the 
uv- Plane Correction 

Before continuing, let's pause to develop an intu- 
ition of what is happening here. Consider a simple ex- 
ponential, 

x = e 71231 w -.53232 - .84654i. (4.1) 

Suppose that we want to approximate this exponential 
using a Taylor expansion, 

.-fe^. (4.2) 

It makes sense that a decent approximation to the origi- 
nal exponential should be possible by truncating this sum 
after a finite number of terms. But how many terms are 
necessary? 
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Fig. 3. — The results from correcting to first order in the uv- 
plane. Top Panel The real space sky (the axes have dimensions 
of radians): After correction to first order, the star remains at its 
perturbed location. The most important feature of this graph is the 
intensity of the sky. In both the pure and corrupted skies of figure 1 
this intensity had been 1.0 at one pixel and 0.0 at all others. Here, 
it takes a maximum value of approximately twice that. Bottom 
Panel The visibility space sky (the axes have dimensions of inverse 
radians): This plot has two important features: 1) The maximum 
in the intensity, which has now increased by an order of magnitude 
near the extreme lower left and upper right corners compared to its 
value in the pure visibility sky (figure 2) and corrupted visibility 
sky. 2) The border around the edge, which is caused by the non- 
local nature of the u-u-plane correction (see equation 2.8) combined 
with the finite nature of our numerical liD-plane. More specifically, 
points in this border strip are not entirely corrected because part 
of their correction requires points in the ™-plane outside of the 
numerical grid. 

Let's first consider the zeroth order approximation, in 
which only the n = term is kept, 

x = 1. (4.3) 

Notice that the zeroth order approximation gives us the 
right magnitude of 1, but all information about the phase 
has been lost. If we instead correct to first order, we 
obtain 

zi = l + 7123i (4.4) 

This is not even close to the correct answer - not only 
does this not contain the correct phase, but the magni- 
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Fig. 4. — The results from correcting to second order in the uv- 
plane. Top Panel The real space sky (the axes have dimensions of 
radians): Notice that the maximum in the intensity is now approx- 
imately 12 times that of the pure and corrupted skies of figure 1. 
(Recall that the real space intensity after the first order correction 
was only 2 times too large). Bottom Panel The visibility space 
sky (the axes have dimensions of inverse radians): Notice that the 
maximum in the visibility sky intensity is 10 times the maximum 
for the first order correction (figure 3, bottom panel) and 100 times 
the maximum for the pure simple visibility space sky (figure 2). 



tude is now not even close to being correct. Correcting 
to second order gets us even further from the correct an- 
swer, 

x 2 = -25368563.5 + 7123i (4.5) 

This pattern continues for higher orders as well. In 
fact, the approximation won't begin to look decent un- 
til ri w 7123. Even more relevant to our observations 
in the previous section, notice that adding subsequent 
terms to the approximation does not necessarily make 
the approximation better until n ks 7123. Before then, 
adding subsequent terms actually makes the approxima- 
tion worse. Drawing from these observations, we expect 
that the trend we've seen so far is the result of under 
correcting in the ww-plane, and that by going to more 
and more orders we will eventually obtain a decent cor- 
rection. 
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Fig. 5. — The results from correcting to fifth order in the uv- 
plane. Top Panel The real space sky (the axes have dimensions of 
radians): The maximum in the intensity continues to increase com- 
pared with the pure and corrupt skies (figure 1) and lower order 
corrections (figures 3 and 4, both top panels), and is now approxi- 
mately 120 times too large. Bottom Panel The visibility space sky 
(the axes have dimensions of inverse radians): The maximum in 
the intensity continues to increase near the lower left and upper 
right extremes, and is 1000 times larger than the maximum in the 
pure visibility sky (figure 2). 

4.4. Higher Order Corrections, 10+ orders 

We now verify that this intuition is correct by 
studying higher order corrections. If correct, we expect 
to see the results gradually get better. We now consider 
the fifteenth order correction, n max = 15, as shown in 
figure 7. The real space sky after 15 orders (top panel) 
is now only 500 times too intense (versus 600 for 10 or- 
ders), while the visibility space sky (bottom panel) is still 
approximately 10 4 times too intense. From this it is un- 
clear that things are getting better, but in the very least 
the intensities are not getting worse. The shape of the 
source, however, continues to grow further from a point 
source. 

Moving on to n max = 20 (figure 8) is a bit more re- 
assuring. The maximum in the real space sky intensity 
is now approximately only 50 times its actual value (top 
panel) , although the maximum in the visibility space in- 
tensity is still four orders too high (bottom panel). The 



Fig. 6. — The results from correcting to tenth order in the uv- 
plane. Top Panel The real space sky (the axes have dimensions 
of radians): The upward trend in the maximum of the intensity 
increases, and the maximum in the intensity is now approximately 
600 times its value in the pure and uncorrected corrupt skies (figure 
1). In addition, the star, initially a single-pixel point source, is now 
spread over an appreciable number of pixels. Bottom Panel The 
visibility space sky (the axes have dimensions of inverse radians): 
Likewise, the maximum in the visibility space intensity continues 
to increase in two of the extreme corners, and is now approximately 
10,000 times what is was in the pure visibility space sky (figure 2) 
and uncorrected corrupt sky. 

gradual improvement continues when we skip ahead to 25 
orders (figure 9). The real space intensity is now only 4 
or 5 times too large (top panel) , while the visibility space 
intensity has now dropped to 1000 times too large (bot- 
tom panel). The shape of the source, however, continues 
to grow worse. 

Skipping ahead next to 30 orders of correction shows a 
dramatic improvement. The top panel of figure 10 shows 
that the single source now appears to be a single source 
of the right order of magnitude in intensity. And in addi- 
tion, the most intense pixel is now located exactly where 
it was for the pure sky, so the uu-plane correction has 
(at least in terms of location of the max) successfully 
corrected for the shift by the ionosphere. The bottom 
panel of figure 10 shows that the maximum in the visi- 
bility space sky intensity (which, as always, occurs near 
the edge of the grid) is now only an order of magnitude 
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Fig. 7. — The results from correcting to fifteenth order in the 
uv-plane. Top Panel The real space sky (the axes have dimensions 
of radians): The maximum in the intensity has decreased relative 
to its 10th order counterpart (figure 6, top panel), from approx- 
imately 600 to approximately 500 times too large relative to the 
pure and uncorrected corrupt skies (figure 1), but the source con- 
tinues to spread in extent. Bottom Panel The visibility space sky 
(the axes have dimensions of inverse radians): While the real space 
sky showed a decrease in intensity compared to the 10th order cor- 
rection, the maximum in the visibility sky intensity is still approx- 
imately 10,000 times too large compared to its pure counterpart 
(figure 2). While the visibility sky intensity does not look much 
better in this sense, it at least does not appear to be getting worse. 
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Fig. 8. — The results from correcting to twentieth order in the 
nu-plane. Top Panel The real space sky (the axes have dimensions 
of radians): The maximum in the intensity continues to decrease 
now as we continue to higher orders of correction. The maximum 
is now only approximately 50 times its value for the pure and un- 
corrected corrupt skies (figure 1), down by a factor of 10 compared 
to the fifteenth order correction (figure 7, top panel). But while 
the maximum intensity looks better, the star has spread even fur- 
ther in extent. Bottom Panel The visibility space sky (the axes 
have dimensions of inverse radians): The maximum in the visibil- 
ity space sky is comparable to the tenth (figure 6, bottom panel) 
and fifteenth (figure 7, bottom panel) order values at approximately 
10,000 times larger than the pure (figure 2) and uncorrected cor- 
rupt visibility space skies. 



too big. It appears as if we've gone over the hump, and 
are now on our way to decent results. 

The correction to 35 orders shows minor improvement 
in real space (figure 11, top panel). In the visibility space 
sky (figure 11, bottom panel), however, the entire grid 
now has the correct order of magnitude of 10 -6 , includ- 
ing the most extreme pixels. We therefore now see some 
of the finer patterns caused by the non-locality of the 
correction and finite nature of the uv grid (as mentioned 
previously in section 4.2), which had previously been hid- 
den by the extreme intensities at the corners. 

It should be noted that the most intense pixel in this 
fully corrected real space sky in the top panel of figure 
11 still lies at the location of the single source in the 
original, pure sky. In other words, the Mf-plane correc- 



tion has successfully shifted the reflected source back to 
its initial position, at the cost of minor spreading over a 
few neighboring pixels. This spreading, which cannot be 
eliminated by correcting to still higher orders, is caused 
by the finite nature of the wii-plane and thus cannot be 
avoided. 

It turns out that corrections to higher orders show neg- 
ligible improvement over the correction to 35 orders, so 
the resulting skies, identical to those of figure 11, are not 
shown. 

So far we have only used a particularly simple sky with 
one star. Figure 12 shows (the absolute value of) the 
residual between a more complicated pure sky with 10 
stars, and the real space sky obtained after perturbing 
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Fig. 9. — The results from correcting to twenty-fifth order in the 
™-plane. Top Panel The real space sky (the axes have dimensions 
of radians): The maximum in the intensity of the real space sky 
continues to decline, and is now only 4 or 5 times larger than the 
actual pure simple sky value (figure 1). The star, however, still docs 
not resemble a point source. Bottom Panel The visibility space sky 
(the axes have dimensions of inverse radians): The maximum in 
the intensity of the visibility space sky is now beginning to decline, 
and is now only 1000 times its value in the pure visibility space 
sky (figure 2), compared to 10,000 times too large for the twentieth 
order correction (figure 8). 



this pure sky and then correcting in the uu-plane to 40 or- 
ders. For comparison's sake the perturbation used here 
was the same ionospheric mode used above to perturb 
the simple sky of figure 1 (ie, that used throughout this 
section) . This figure shows the kinds of errors we may ex- 
pect from the uw-plane correction. The residual from the 
star on the top right shows a light cross pattern, indica- 
tive of a spreading of the source caused by the process of 
perturbing the star and then applying the uv-plane cor- 
rection. Places in the plot with two consecutive pixels of 
high intensity represent stars which were not shifted back 
to exactly the same pixel that they started at, but rather 
to a neighboring pixel. Recall that the stars are typically 
initially shifted 3 or 4 pixels by the ionosphere, so the uv- 
plane correction is still providing some improvement with 
these stars. 

5. ANALYTIC ESTIMATE FOR N MAX 
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Fig. 10. — The results from correcting to thirtieth order in the uv- 
plane. Top Panel The real space sky (the axes have dimensions of 
radians): The single source from our simple sky (figure 1) now once 
again looks like a simple source, but with its intensity spread over 
a few neighboring pixels. The maximum in the intensity is now of 
the same order of magnitude as that of the pure sky (figure 1), and 
is located at the same pixel. Bottom Panel The visibility space sky 
(the axes have dimensions of inverse radians): The maximum in 
the visibility space intensity, which is located at the lower left and 
upper right corners as always, is now only an order of magnitude 
higher than that of the pure visibility sky (figure 2). 

There are two main points to be taken from our 
analysis so far. First, from section 3 we learned that ei- 
ther needing too many modes to model the ionospheric 
correction or too many orders of correction leads to an 
unreasonable number of numerical calculations. Second, 
from section 4 we learned that under correcting in the 
uu-plane is a huge mistake and a lot worse than not cor- 
recting at all. Hence the dilemma: choosing n max too 
small leads to the destruction of the data, while choos- 
ing n max too large leads to a computationally infeasiblc 
problem. It is therefore advantageous to develop a the- 
oretical prediction of how many orders of correction are 
necessary. As it turns out, the result will lead to a few 
tricks which make the problem more reasonable. 

5.1. Finding an Upper Bound on the Error U(u), and 

Umax 
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Fig. 11. — The results from correcting to thirty- fifth order in the 
MD-plane. Top Panel The real space sky (the axes have dimensions 
of radians): This real space sky shows minor improvements over 
the real space sky corrected to 30 orders in the MD-plane (figure 
10). For both, the source has been successfully shifted back to 
its correct starting pixel, as shown in the top panel of figure 1, 
at the cost of minor spreading to the neighboring pixels. Bottom 
Panel The visibility space sky (the axes have dimensions of inverse 
radians): The entire visibility space sky (including the extremes) 
now has an intensity of the order 10 — 6 , which is the order of the 
intensity for the pure visibility space sky (figure 2). Notice that 
since the most extreme pixels are no longer much more intense than 
the others, we now see the pattern caused by the non-locality of 
the correction and finite nature of the grid (see section 4.2). 

If we correct to only n max orders, then the magni- 
tude error E{u) in our result must be the absolute value 
of the sum of all the terms we left out; more specifically, 

oo / / 2M / ^ j \l \ 

E ^= E E n^r 1 - W+EU) 

(5-1) 

The steps leading to an upper bound U(u) on this error 
may be found in the appendix, section B. The result is 



U(u) = 
where 



ImaxV(u) 



n (n + 2M - 1)! 



v(u) = max { \c q u- d q \ } any q, 



^rr > E(u), 

n!(2M-l)! (( 5 2 ? )!) 2M ~ 

(5.2) 
(5.3) 




Fig. 12. — 10 Star Sky This plot shows the residual between a 
pure real space sky with 10 stars, and the real space sky obtained 
after perturbing this pure sky and then correcting in the MD-plane 
to 40 orders. The ionospheric perturbation was the same used for 
all the other figures in this section, figures 1 to 11. This figure shows 
the errors typical with the iro-plane correction. The top right cross 
represents a star shifted back to the correct position, but with an 
intensity spread over a few pixels. Stars with two consecutive high 
intensity pixels represent stars shifted back to pixels adjacent to 
their correct starting pixels. (Typically, stars are shifted 3 or 4 
pixels by the ionosphere, so the u?;-plane correction has provided 
some improvement with these stars as well). 



max { |/(u) } any u, 



(5.4) 



and M is the number of modes, as always. [See equa- 
tions^. 9 and 2.10 for reminders on how^effective modes 
(c q , d q ) are related to actual modes (a m , b m )]. This is our 
final result for a strict upper bound on the error. Unfor- 
tunately, this formula is not too enlightening. In order 
to obtain a theoretical estimate for n max , we must make 
a few further approximations. 

As one would expect, the optimal value of n max , which 
represents the number of orders necessary to obtain some 
level of accuracy in the ww-plane, is dependent upon the 
level of accuracy desired. To quantify this, define I nmax 
to be the value of the intensity in the uu-plane after cor- 
recting up through n = n max , and I ac tuai to be the value 
of the intensity in the uw-plane that one would obtain by 
employing the full correction and not truncating the sum 
(ie, n max — oo). (It should be noted that I ac tuai here also 
assumes a wv-plane infinite in extent. This will have ef- 
fects seen later.) The fractional error / in the uu-plane 
correction caused by truncating the sum is then 



/ = 



|-^i maa . -^-actual \ 
| ^actual | 



(5.5) 



In the appendix (section C) you will find the steps leading 
up to a theoretical prediction of the value of n max at a 
given u in the itv-plane necessary to obtain a fractional 
error less than or equal to P if given the ionospheric 

effective modes distorting the sky, (c q ,d q ). The result is 
that the optimal value of n max is estimated by 



min { n } such that Pnl > (2Mv(u)) r 



(5.6) 
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Fig. 13. — Theoretical Predictions for n m ax- Top Panel The 
theoretical predictions for n max shown graphically as a function of 
u for the simple sky explored in section 4 (and shown in figure 1). 
This plot has P = .1 (see equations 5.5 and 5.6) . The axes are 
measured in inverse radians, and the third dimension with color 
represents the calculated value of n ma x ■ The preferred direction is 
defined by the direction of the ionospheric perturbation (tt/4 with 
respect to the u x axis). Notice that the value of n max increases 
as u moves radially outward along this direction. Bottom Panel 
The residual between the theoretical prediction shown in the top 
panel and the numerical calculations of n max for the same sky. 
For a bulk of the MD-plane, the prediction is quite good (within 
5). The cool colored areas represent pixels that never obtain the 
desired fractional error due to the non-locality of the correction 
(see section 4.2). 



where M is the number of modes and v(u) is defined as it 
was above in equation 5.3. This formula is a little hard 
to digest, so some values for n max given Z = 2Mv(u) 
and P are provided in table 5.1. For the example sky and 
perturbation used in section 4, the theoretical predictions 
for the number of orders necessary is shown graphically 
in the top panel of figure 13. This figure has P = .1, 
although n max does not change too significantly when 
varying P, as is seen in table 5.1. Recall that the single 
distorting mode d q is in the 7r/4 direction, which defines 
the favored direction seen in this figure. 

As a check of these theoretical predictions, we used 
MATLAB to numerically compute the number of orders 
necessary to obtain the desired fractional error of P = .1. 
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The bottom panel of figure 13 shows the difference be- 
tween these computational results and the theoretical 
predictions shown in the top panel of figure 13. More 
specifically, it represents the number of orders of cor- 
rection theoretically predicted minus the number found 
computationally. This figure suggests that for a bulk 
of the uv-plane, the theoretical prediction is quite accu- 
rate, predicting the number of orders to within 5. Near 
the extremes, however, the finite nature of the uu-plane 
causes problems (remember that the theoretical estimate 
assumed an infinite ww-plane). In fact, the cool-colored 
pixels near the corners are pixels which never obtained 
a fractional error of P < .1 (The numerical code cutoff 
after 50 orders; all points with fractional errors too high 
at that point were assigned a value of 50 orders). 

The theory predicts that about 35 orders are re- 
quired to correct at the most extreme points in the uv- 
plane, which is what our previous numerical computa- 
tions found. An important feature of the theoretical pre- 
dictions shown in figure 13 that is characteristic of all 
skies is that the necessary number of orders of correction 
varies with u, and in particular it increases as u increases 
along the direction of the mode. Therefore, points closer 
to the origin are corrected in less orders than those fur- 
ther away. 

5.2. The Strongest Mode and Significant Modes 
Approximations 

The accuracy of the theoretical prediction here is 
in no small part due to the existence of only one iono- 
spheric mode in our simple sky model. This reason for 
this is that the above theoretical estimate (equation 5.6) 
was derived from an expression for the upperbound on 
the error (equation 5.2) which assumes that all the iono- 
spheric modes in the sky are as strong as the strongest 
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mode at u and add constructively (which may be seen in 
the appendix, section B, near equation B.16). As such, 
the result is not a bad prediction for only one distorting 
mode, but tends to (perhaps significantly) overestimate 
the necessary number of orders for multiple distorting 
modes. In short, the above mentioned theoretical esti- 
mate may perhaps be more accurately called a theoret- 
ical overestimate. Given the results seen in section 4 
(more specifically, the terrible consequences of undercor- 
recting in the iw-plane), this was done intentionally to 
ensure that our uw-plane was adequately corrected. Still, 
it may be useful to obtain a more accurate estimate of 
the number of orders necessary. 

One such estimate would be a strongest mode approx- 
imation, in which we assume that at any given u, the 
only significant contribution comes from the strongest 
mode at that point. The contributions from the other 
modes are assumed to be weak and negligible. This ap- 
proximation ultimately boils down to setting M = 1 in 
the final equation determining n max from the previous 
section (equation 5.6). This approximation may provide 
a more accurate estimate of n ma x, but it also runs a 
high risk of underestimating the correct number of or- 
ders, which should be avoided if possible. An alternate 
approximation would be a significant modes approxima- 
tion, in which only modes at a given u with strengths 
within a certain critical fraction of the strongest mode's 
are included in the value for M used in equation 5.6. 
With a closer study of perturbations from more realistic 
ionospheric modes, it may be possible to set this critical 
fraction in such a way as to fairly accurately predict the 
number of orders necessary. 

6. TWO METHODS FOR MAKING THE (7V-PLANE 
CORRECTION MORE FEASIBLE 

The above analysis suggests two methods for mak- 
ing the ww-plane correction less time intensive: 1) The 
points furthest out in the uw-plane take the most time to 
correct. Eliminating them decreases computation time, 
but at the cost of resolution in the real space sky. 2) Dif- 
ferent points in the uu-plane require different numbers 
of orders of correction, so write a code that corrects to 
different numbers of orders depending on the point in the 
uu-plane. (In other words, don't waste time correcting 
to 35 orders near the origin when 2 is enough.) 

6.1. Method 1: Eliminating the Extremes of the 
uv -plane 

For this method, we eliminate the problems caused 
by under correcting at the extremes in the ww-plane by 
setting the values at those extremes to 0. Take, as a vi- 
sual example, figure 14, in which we have set the values of 
the pixels in the 25 diagonal rows from the corner to zero 
(we shaved 25 pixels from the corner) . As a reminder of 
the real space sky after 25 orders without edge shaving, 
consider the top panel of figure 15 (a reproduction of 
figure 9, top panel). Notice that the intensity is approxi- 
mately 4 or 5 times too high at the brightest points and, 
even worse, our single point source has turned into some 
sort of supernova explosion. Compare this to the edge- 
shaved version of the real space sky (figure 15, bottom 
panel), in which we see a sky that looks almost identi- 
cal to our fully corrected sky after 35 orders (top panel, 
figure 11). To see this more clearly, consider figure 16, 



Visibility Space Sky, C: 25 orders, S: 25 pixels, corner 




Fig. 14. — The Edge Shaving Technique This is a visual exam- 
ple of the edge shaving techinique, in which the pixels in the 25 
diagonal rows from the upper right and lower left corners (ie, the 
pixels which require the most orders of correction) have had their 
values set to zero (or, are shaved). Notice that eliminating these 
ultra intense pixels reveals some of the finer patterns caused by 
the u^-plane correction, as also seen in the correction to 35 orders 
(figure 11, bottom panel). Compare this to the uv-plane created 
by correcting to 25 orders and not edge shaving (figure 9, bottom 
panel), in which the intense pixels from the corner dominate the 
others in the nu-plane. 

which represents (the absolute value of) the residual be- 
tween the real space sky corrected to 40 orders with no 
edge shaving and the real space sky corrected to 25 orders 
with edge shaving. From this figure we see that the result 
of the edge shaving was to create a small spread around 
the star, but of an intensity about an order of magnitude 
lower than the maximum intensity of the star. 

The relative success of this scheme leads to the ques- 
tion of how low we may push the number of correction 
orders when edge shaving is introduced. Figure 17 shows 
the result of only correcting to 15 orders, but shaving 75 
rows of pixels from the corners in the uu-plane. Without 
edge shaving, the real space sky corrected to 15 orders 
had a maximum intensity of about 500 (figure 7, top 
panel). Now, the total intensity is approximately 1 as it 
should be, but it is spread over a number of pixels. So 
while the results are a dramatic improvement over what 
they had been, for the sake of precision it might be a 
good idea to correct to higher orders and shave less. The 
moral: The process works, but be careful about trying 
to shave too much. 

6.2. Method 2: Correcting to Different Orders at 
Different u 

For this method, we correct to different numbers of 
orders at different u. To test this method, we re-wrote 
our MATLAB code so that the number of orders of cor- 
rection at a given u was determined by the theoretical 
estimate from section 5 (More specifically, equation 5.6). 
We then corrected the distortion for the same simple sky 
used in section 4. Figure 18 shows the residual between 
the real space sky corrected to 40 orders at all points in 
the itu-plane and the real space sky corrected to different 
orders in the ui>-plane. Not surprisingly, the residual is 
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Fig. 15. — The results of the edge shaving technique. Top Panel 
A reproduction of the top panel of figure 9, shown as a reminder 
of the real space sky that results by correcting to 25 orders in the 
«D-plane and not edge shaving. Bottom Panel The real space sky 
that results from correcting to 25 orders and shaving 25 rows of 
pixels from the corner in visibility space (ie, the real space sky cre- 
ated by inverse Fourier transforming the tre-plane shown in figure 
14). Notice the drastic improvement: The sky now looks like a 
single source of approximately the right intensity and at the right 
location, but with a small spreading of intensity to neighboring 
pixels. 



incredibly small- 5 or 6 orders of magnitude less than 
the maximum in the intensity of the source. (This is, 
of course, another sign that the theoretical prediction of 
the number of orders of correction is pretty good). How- 
ever, this new MATLAB code presented a small problem: 
MATLAB is so much better at manipulating matrices 
than running for-loops that this second code, which the- 
oretically requires less computation, takes approximately 
15 times as long to run. Of course, if the Mf-plane correc- 
tion is eventually used in MWA, a programming language 
more adept at loops will undoubtedly be used, and this 
method will potentially save time. 

7. AN ALTERNATE CHOICE FOR 59 
Previously, we have assumed that 59 took the form 




Fig. 16. — This is the residual between the real space sky created 
by fully correcting to 40 orders with no edge shaving and the real 
space sky created by fully correcting to 25 orders and then edge 
shaving 25 rows of pixels (shown in the bottom panel of figure 
15). Notice that the two skies seem to match pretty well, with 
with maximum disparities an order of magnitude lower than the 
maximum intensity of the source. 




Fig. 17. — The real space sky that results from correcting to 
only 15 orders, and then shaving 75 diagonal rows of pixels in the 
uv-pl&nc. Comparison of this real space sky to that obtained by 
correcting to 15 orders and not edge shaving (see the top panel of 
figure 7) shows remarkable improvement. The maximum intensity 
has decreased by a factor of 1000, and is now of the right order of 
magnitude. On the down side, however, the source is now spread 
among a fair number of pixels, and might not posses the precision 
one desires. 



A I 



59' = Re{J2 
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(7.1) 
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Given the relative complexity of the results above, it pays 

to investigate an alternate choice of 59 . The main rea- 
son that this was chosen was because it allowed for an 
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t x =J2 kl k, (7.7) 

k 

and the sums over l m and Ik are restricted so that 
J2k l k = n x and J2 m l m = n y . 
While we've confined the analytic solution of this to 

the appendix, it should be mentioned that en route to 

-> i 

this solution the exponential in 59 from equation 7.2 
was Taylor-expanded (which, you may recall, was also 

the case for the other choice of 59 as a sum over sinu- 
soidal modes, and is explicitly shown in equations 2.1 
and 2.2). In other words, this solution is likewise char- 
acterized by the double expansion mentioned previously 

at the end of section 2: one expansion over (n x , n y ) from 

— i 

expanding the exponential in 59 , and one expansion in 
(^fc)'m) resulting from our model for the ionosphere. As 
such, this solution shows many of the unfortunate char- 
acteristics of the sinusoidal choice. In particular, there 
are still restricted sums which contain a number of terms 
comparable to that calculated in section 3, so this choice 
has the same problem of making the uw-plane correction 
unreasonable if too many orders of correction or iono- 
spheric modes are needed. 

In addition, the correction is still not strictly local: 
numerical computation of derivatives requires neighbor- 
ing pixels, with higher orders requiring more neighbors. 
Moreover, numerical computation of derivatives for a fi- 
nite data set introduces its own set of additional errors, 
and thus makes this choice much less appealing than the 
previous one with sinusoidal modes. 

The ugliness of both of these choices ultimately stems 
from the inability to solve for I(u) analytically with- 

out expanding the exponential containing 59 . Unless 
a choice is found which may be solved analytically with- 
out this first expansion, it is doubtful that a better choice 
than the sinusoidal modes will be found. 

8. A MORE REALISTIC IONOSPHERE 

The ionospheric distortion presented throughout 
this paper was chosen because its particularly strong 
nature accentuated the subtleties of the ww-plane cor- 
rection. The strength of this mode made the uu-plane 
correction appear computationally infeasible: any iono- 
sphere which required even 10 of these modes to accu- 
rately model would require too much computation (see 
section 3). However, an ionospheric mode which shifts 
sources on the sky by tens of arcmin is somewhat un- 
realistic. We conclude by considering a more realistic 
distortion. 

The computational feasibility of the uw-plane correc- 
tion is determined by the largest value of n max (given by 
equation 5.6) required for any u. To calculate this for a 
realistic sky, we need to know the largest possible value 
of v(u). It is possible to cast the largest value of v(u), 
which we label v max , in a form which better elucidates 
its physical significance. More specifically, notice that 

v(u)= max { \c q u ■ d q \ } any q (8-1) 
= max { |cg|u||dg| cos(ot)| } any q, (8.2) 

where a is the angle between u and d q . For an arbitrary 
choice of u, this cosine term may be significant. However, 



Fig. 18. — The residual between the real space sky created by 
correcting to 40 orders at all points in the tt-u-plane and the real 
space created by correcting to varying orders in the iii)-plane, as 
determined by the theoretical prediction for n max (equation 5.6; 
shown graphically for the simple sky of section 4, which was also 
used here, in the top panel of figure 13). Notice that the maximum 
difference is 5 or 6 orders of magnitude less than the maximum in 
the intensity (Recall that the correction to 40 orders is essentially 
the same as that to 35 orders, which is shown in the top panel of 
figure 11). This suggests once again that the theoretical prediction 
for the number of orders is fairly accurate. 



analytic solution for I(u), where 

m\= j f^L I (fe'^'-^-^e-^'im, (7.2) 

J (2tt)^ J 

as shown previously in equation 1.17. Without an an- 
alytic solution to these integrals for I(u), the uu-plane 
correction becomes more computationally intensive than 
inverse Fourier transforming to real space, correcting for 
the ionosphere there, and then Fourier transforming back 
to visibility space. Therefore, an analytic solution is re- 

quired for any choice suitable choice of 59 . Unfortu- 

nately, choices for 59 which allow such analytic solutions 
are hard to find. There is, however, at least one other 
such distortion: a polynomial expansion, given by 

K M 

59' = Y / ake' x k +Y, b ™ 9 T> ( 7 - 3 ) 

k— 1 m— 1 

where K and M are the number of terms in the 9 X and 
9 y directions, respectively, necessary to accurately model 
the distortion by the ionosphere. The au and b m shown 
here are real. The analytic derivation of I(u) is given in 
the appendix, section D. The result is 

oo oo / / 

HQ = E E E E -(-i««)'*--*-+ l (-<u v ) n »- t » +1 (7.4) 

"« n v lm h 

m=l lm - h=l lk ' du * 9u V 

where 

t y=y2 ml m, (7.6) 
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if we wish to calculate v max , we may set a — and \u\ 
= Umax, where u max is the greatest distance from the 
origin in the uw-plane that our antenna's uw-coverage 
allows. This is valid for MWA because the wv-coverage is 
approximately circular, so that the strongest ionospheric 
distorting mode is guaranteed to lie along a direction 
which possesses this maximum displacement in the uv- 
plane. With these changes, and based upon our previous 
definitions of the effective modes {c q , d q } in terms of the 

actual modes {a m ,b m } (see equations 2.9 and 2.10), we 
may write v max as 

v max = max { |— a m bm\Umax } any m. (8.3) 

But, 

\a m b m \ = \ia m b m e- lb ™- '\ (8.4) 
> \Re{iaJ m e-^ S ')\ = \66'\ (8.5) 

Let's define 

| = max { | | } any m. (8.6) 

In words, \59 max \ is the maximum deflection caused by 
a single mode that we might observe. In addition, let B 
denote the length of our antenna array's longest baseline. 
The maximum wu-planc displacement u max is then this 
length measured in units of the wavelength A that our 
antenna is detecting, 

Umax = — (8.7) 

[The extra factor of 2ir is the result of our convention for 
Fourier transforms (see equation 1.7), which differs from 
that conventionally used in radio astronomy] . With these 
substitutions, 

1 2,ttB 

Umax = - \S0 max \ - . (8-8) 

Therefore, in terms of these parameters, the number of 
orders of correction necessary is (adapted from equation 
5.6) 

n m ax = min { n } such that Pn\ > (M\89 max \^^-) n , 

(8.9) 

where P is (as before) the fractional error desired for 
the correction. It should be noted that this form is only 
valid for determining the largest n max among all u. For 
calculating n max for a particular u, equation 5.6 must be 
used. 

For MWA, a typical frequency detected will be about 
140 MHz, corresponding to A « 2.143 meters (this rep- 
resents the 21 cm emission for a red shift of z « 9.2.) 
We expect the ionosphere to deflect such a wave approx- 
imately 0.6 arcmin = 1.75 * 10~ 4 radians ((3), value is 
for the night). If we consider baselines of approximately 
400 meters, then 

Z full ee M\66 max \^ w .205M. (8.10) 

If given the number of modes necessary to accuately 
model the ionosphere M (which is as of yet undeter- 
mined), then n max for the full array may be determined 



from table 5.1 by substituting Zf u u for Z. As an ex- 
ample, if a fractional error of P = .1 is desired and 
M = 20, then Zf u u ss 4 and table 5.1 shows that 12 
orders of correction are necessary. Whether such a re- 
sult is computationally feasible is dependent upon how 
much time is alloted for the ionospheric correction and 
the quality of the computers used. As such, it may not 
be determined here. What is clear, however, is that such 
a correction is not obviously ruled out on computational 
grounds (especially if a technique such as edge shaving 
is used to reduce n max ). [Quick aside: Edge shaving al- 
ters the above results by substituting the largest |u| left 
unshaved in place of u max in the above calculations.] It 
should be noted that the values that went into calcu- 
lating Zf u u above were estimates, and certainly not set 
in stone. In particular, we once again emphasise that 
throughout this paper we have remained ignorant of the 
details involved in the expansion of the ionosphere, and 
have no knowledge of how many modes M are necessary 
to sufficiently model the effect of the ionosphere. In addi- 
tion, the ionospheric deflection l^ma^l is proportional to 
A 2 , with longer wavelengths experiencing greater shifts 
(2). Therefore, we expect longer wavelengths than the 
above to require more orders of correction and shorter 
wavelengths, fewer. If Zf u u is in fact lower by a factor 
of 5, for example, then n max « 4 and the wv-plane cor- 
rection is certainly a viable candidate for correcting the 
ionosphere. In particular, if the strongest mode approx- 
imation discussed in section 5.2 turns out to be a good 
approximation, then even with baselines of 1.5 km we 
may expect a good correction after only 4 orders for the 
wavelength given above. On the other hand, if Zf u u is 
raised by a factor of 5, then n max « 54 and the uu-plane 
correction is clearly computationally infeasible for any 
reasonable value of M . 

9. POTENTIAL COMPUTATION SAVER: UPDATING THE 
A T MATRIX 

As stated previously, correcting for the ionosphere 
in the uw-plane entails multiplying the perturbed data 
I(u') by the ionospheric correction operator A T (u,u'), 

I(u;t) = A T (u,u';t)I(u';t), (9.1) 

where we have now made the time dependence of these 
quantities explicit. The above may be thought of as a 

matrix equation, where I{u';t) and I(u;t) represent our 
uncorrected and corrected (respectively) data arrays, and 
A T (u,u';t) represents a correction matrix. The entries 
of this correction matrix are calculated by the appropri- 
ate binning of the coefficients in our correction equation 
(reproduced from 2.8), 

'(**)= E E IT {CqU l i - U(a+Ew)> 

n = ll,l2,---,hM \4=1 " / 1 

(9.2) 

where the "coefficients" are the quantities preceeding 
I(u + ^2 q lqdq;t) on the right hand side of the equa- 
tion. If the timescale within which one wishes to recal- 
culate the effect of the ionosphere is small compared to 
the timescale within which the ionosphere significantly 
changes, then it is possible that the correction matrix 
A T (u,u';t) has changed very little from that previously 
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calculated. More specifically, if one wishes to calculate 
the correction matrix at a time ti shortly after having 
calculated it at time to (ie, if \t\ — t \ « r where r is 
the time scale of significant change in the ionosphere), 
then 



A T (u,w';ti) w A T (u,u ;t ) + 5A T (u,u;t ) (9.3) 



where i5A T represents a small correction matrix. In this 
regime, it is computationally much easier to calculate the 
small correction <5A T and add it to the previously cal- 
culated A T (u, u'\ to) then to calculate A T (u, u'\ t\) from 
scratch. Therefore, in such a scenario the process of up- 
dating the correction matrix is computationally favor- 
able. 

The first order correction 5A T may be analytically cal- 
culated as follows: Let's assume that^we model the iono- 
sphere using the effective modes c q , d q that we have been 
throughout this paper (see equations 2.9 and 2.10), and 

that the d q represent part of a fixed Fourier basis while 
the c q are our fitting parameters. (In other words, the 

d q are fixed and time independent, while the c q fluctuate 
with time). Let's assume that we've calculated the cor- 
rection at a time to- More specifically, assume that for 
all u we've calculated and stored all the relevent terms 
in the correction equation for I(u; to), 



e e (n ^ (t f, rf > )( 9 .4) 

n=0 Ii,12, — Mm \<?=1 / 

x /(u + EM?;*o), (9.5) 



where the time dependence of /, /, and c q is now explicit. 
Now we want to correct for the ionosphere at a later time 
t\. To compose the correction matrix at time ti, we could 
start with the full correction formula, 



'1M 



e e n 



(Cq(tl)u ■ d g y° \ 

n=0 h,l 2 ,...,l2M \<J=1 J 

x j(t?+Ei«4;*i). ( 9 - 7 ) 



and then construct the new matrix A T (u, u'; t\). In- 
stead, however, let's assume that we're in the regime of 
small ionospheric changes, so that to first order 



c q (ti) w c q (t ) + Ac, 



(9.8) 



rection and only keeping terms to first order we obtain 



' «c q (t ) + Ac q )u ■ d q ) l « 



i{u M )= E E n 

n=0 h,l2,---,hM \9=1 

x /(« + E^4;*i) 



u 



n=0 h,h,---,hM \9=1 

X /(u + EM<,;ti) 
q 

e e n (Cg(f0) M q } i+ e 



n=0 l 1: l 2 ,---,l2M \9=1 

x /(« + E^4;*i) 



1=0 h,h,---,hM \9=l 9 



E e (n^^i^Ev; (1 )- 



n mox / /2M , ,. r w \ / 2M . 



e e n 



where Ac g is small. Substituting this into the full cor- 



= h,l2,-,hM \q=l ) \q=l 

X I(u + y^J q d q ;ti). 

q 

In this final equation, the coefficients in this first group 
of terms exactly replicate those coefficients from time 
to (sec equation 9.4). In other words, these terms 
represent the previously determined correction matrix 
A^(u, u'; to)- The second group represents the small ad- 
justment 5A T to the correction matrix at to- Notice 
that these terms have been there own maximum cutoff 
for ra, labeled as n' max in the above equation. If the 
quantity Ac q /c q (to) is small, as assumed, then the in- 
dividual terms in this second sum are also small, and 
thus a smaller value of n' max is necessary to obtain a de- 
sired fractional error for the intensity. In this case, it is 
computationally favorable to update the correction ma- 
trix rather than derive it from scratch. It should also 
be noted that if n' max < n m ax, then the nonzero en- 
tries of the matrix 6A T form a subset of the nonzero 
entries of the matrix A T (u, u'; to), and thus their sum 
(which represents A T (u, u'; ti)) is equally as sparse as 
A T (u, u'; to). In other words, the matrix A T (u,u';t) 
does not become less sparse through this process of cor- 
rection (a fact which is important for large numerical 
matrix manipulations). 

It is worth mentioning that whether updating the cor- 
rection matrix is a viable method depends on the time 
scales of the changing ionosphere. More specifically, the 
above was calculated keeping only terms to first order in 
Acq/c q (to)- It is possible that for the time scales con- 
sidered higher terms are also necessary, or that Ac q is 
not small compared to c 9 (to); the former situation com- 
plicates the math but does not necessarily outrulc this 
method, while the latter pretty much requires that the 
correction matrix be built from scratch every time. 

10. CONCLUSIONS 

The uw-plane correction only makes computational 
sense if the model for the ionospheric perturbation allows 
for an analytic solution to I(u) (equation 1.17). One such 
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model is a sum over sinusoidal modes (equation 2.3). By 
running numerical codes with this choice, the most im- 
portant result discovered was that under correcting in 
the uv-plane is worse than not correcting at all (section 
4). But in addition to this, correcting to too many orders 
or requiring too many modes to model the effect of the 
ionosphere may lead to a computationally unreasonable 
problem (section 3). To help avoid this issue, a theoret- 
ical estimate of the number of orders of correction nec- 
essary (which agrees well with the sample sky provided 
in this paper) may be used (section 5). This estimate 
reveals that the number of orders of correction necessary 
varies in the ww-plane. This, however, suggests two meth- 
ods for alleviating the problem: eliminating those points 
in the wu-plane which are particularly troublesome at the 
cost of precision for the real space sky (section 6.1) and 
correcting to different orders at different points in the 
■uv-plane (section 6.2). Both techniques prove successful 
and make the problem of correcting in the uu-plane more 
feasible. In addition, depending on how often the iono- 
sphere's effect is updated compared to the timescales of 
change in the ionosphere, it may be compuationally fa- 
vorable to update the previously determined effect of the 
ionosphere rather than rederive its full effect from scratch 
each time. 
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APPENDIX 

The purpose of this appendix is to rigorously derive 
some of the mathematical formulas merely stated within 
the main text. It is included for completeness and for 
the curious reader; no new results are derived. 

A. SOLVING FOR 7(17) FOR A SINUSOIDAL REFLECTION 

In this section, we solve for the unperturbed intensity 
1(H), 

(A.l) 

using a sum over sinusoidal modes for our ionospheric 
deflection 59' , 



M 



59' = Re(J2ia m b m e- lb ^ '). 



(A.2) 



With a bit of algebra (and keeping in mind that b m and 
9' are real but a m is complex), 59' may be written as 

M 



m—1 



■ a mbm ib m -0> _ -KrSra jb„ 

2 2 



59'=J2\ ^e" 6m '"' - ^e lS ™' e ' ) , (A.3) 



where a* m is the complex conjugate of a m . With this 
choice of 59', the expression for I(u) becomes 



J (2tt) 2 J 2n! 



(A.4) 



X ( jr a m u ■ b m e- tim "' - a' m u ■ b m e Sm () ' ) n I(u%A.5) 



Before proceeding, it is convenient to convert the sum- 
mation over m as follows: 

^ a mU- b m c -j$ m .0> _ a* m u-b m c ib m e'^j ^ 

(A.7) 



m—1 



2M 



9=1 



where 



and 



c q = ? 



\a q if q < M + 1; 



if«>M + l. 



(A.8) 



d - I K H q < M + 1; 
9 ~\-b q -M ifg>M+l. (A ' 9) 

With this form, we see that although there are M modes 
distorting the sky, there are 2M terms in the sum. This 

extra factor of 2 comes from the above constraint that 59J 
be real. We will refer to these modes labeled by [c q ,d q 
) as effective modes. Writing the intensity I(u) in terms 
of effective modes gives 

I w=lwr{J** e **~* (A - 10) 



/ 1 \l 



n=0 \q=l 



/(•a'XA.ll) 



The individual terms inside the summation over n may 
be manipulated using the multinomial expansion to give 



/ 2M 



- [^c q u-d q e- 



\q=l 



(A.12) 



= 1 V - 



ll,l2,---,l2M <?— 1 



n (c q u-d q e- iS "- ff y\.ii) 



e-^-E,^', Yl f {CqU ', , 1 ) , (A.14) 



h,h,---,hM 9=1 V 

where E' denotes a restricted sum such that E 9 l q = n 
and l g > 0. The expression for I(u) now becomes 



a) 



(A.15) 



oo / _ 2M , -, r \l a \ 

x E E e- ie ' ^ Sqlq n qU : , imiz) 

n=0i 1 ,i 2 ,...,i 2M « = 1 «■ j 

Because of the above performed multinomial expansion, 

the integral over 9' (once brought inside the summation) 
now takes on the familiar form of a delta-function, and 
is easily performed to yield 



I(u) = j d 2 u' e E <5(«'-s-E <*"«'«) II 



(c q u ■ d q yi 



/(„')• 



(A.17) 

The integral over u 1 is now a simple delta-function inte- 
gral, and its integration gives 

oc / / 2M , - j \l \ 

= £ E II J(g+EW(A-18) 

n=0h,l 2 ,...,l 2M \q=l q ' J 1 
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B. FINDING AN UPPER BOUND ON THE ERROR 

In this section, we estimate the error in the uii-plane 
accumulated by truncating the infinite sum over n after 
the term n = n max + 1- Recall that the total correction 
term is 



e n 

n=0 i1.l2.---.l2M \9=1 



(CgU-dgf" 



l q \ 



I(u + Y,d q l q ). 

(B.l) 



When correcting through n = n max , the magnitude error 
E(u) in the correction is equal to the absolute value of 
the sum over all terms left out. More specifically, 



E(u) 



e e in 

« = 1max + l ll,h,---,hM 



(c q u-d q ) l i 



q i B - 2 ) 

We now attempt to determine an upper bound U(u) on 
this error. To begin, we bring the absolute value inside 
the sum, so that all terms now add constructively, 



2M , _ >; 

>: >: ( i r r '^; ,: )'--••>: 'v.- 

n=n ma ,x + lh,h,---:hM \9=1 q ' J 1 

We expect the magnitude of I(u) will be approximately 
the same at all points in the uu-plane. Denoting the 
maximum value of \I(u) \ for the uncorrected sky as I ma x, 
we find 



E E Imax 

n=n max + l li,l2,---,l2M 



2M 

n 

9=1 



{CqU-dq) 1 " 



l q \ 



> E(i 



(B.4) 

Next, consider the terms in the sum of the form c q u- d q . 
Define v{u) to be the maximum value of \c q u ■ d q \ for 
a given value of u and the effective modes in question. 
Then, 



i1.l2.---.l2M \9— 1 



t, I fi^r H«->". t in') ••• 



I1.I2.---.I2M 
2M 



2M 
[ 

1 ' 



where the equality in the equation above occurs because 
the restricted sum over Z g requires that ^l q ~ n. With 
this substitution, the upper bound on our error function 
becomes 

oo / /2M 1 \ 

E E Iln 

1=1m«i+l il,(2,...,i2M \9=1 *'/ 

(B.7) 

Now concentrate on the inner summation. Define G ac- 
cording to 



G = max 



I assert that 



(2M \ 2M 
Htt where E^ = n - ( R8 ) 
,=1*97 9=1 



G 



,2M • 



(B.9) 



[Aside: Before providing the proof of this, we should 
point out that the fraction n/(2M) is not guaranteed 



to be an integer, and therefore this factorial and the 
ones given hereafter should be taken to be given by the 
Gamma-function T(z), 



z\=T(z + l) 



-f 

Jo 



t z e- l dt\ 



(B.10) 



The proof is quite short: 

I) Start with the given form of G, corresponding to 
l q = n/(2M), for any q. 

II) Any value of JJ f-j which corresponds to an alterna- 
tive choice for the l q under the constraint that ^2 l q = n 
may be obtained by multiplying this value for G by a fi- 
nite number of factors whose magnitudes are all less than 
1. Therefore, G is the maximum. QED. 

An example may be quite useful here: Consider the 
scenario with n = 12 and M = 2. In this case, we assert 
that 

1 



G 



3!3!3!3!' 



(B.11) 



which corresponds to the choice l q = 3 for all four l q . 
Now, let's pick an alternative choice for the l q ; let's say 
h = 5, h = 2, 13 — 4, 1 4 — 1. For this choice we obtain, 



G' 



1 



5!2!4!1!' 
But this may be re-written as 



(B.12) 



G' 



1 



3*2 



3!3!3!3! V4 / \ 4*5 



G * ( factors < 1), 



(B.13) 

where the first factor transforms (^2^3) from (3,3) to 
(2,4) and the second transforms (Zi, Z4) from (3,3) to 
(5,1). Therefore G' < G. 

Our new knowledge of G, when combined with our 
previous determination of the number of terms in the 
restricted sum over l q (see section 3), leads us to conclude 
that 



(# terms) * max | TJ — - 

U=i l «- 



(n + 2M - 1)! 



„!(2M-1)! ( (aifr) |) 



e n 



Il.l2.---.I2M \l = 1 lq ' 



(B.15) 



Therefore, our new upper bound U (u) on the error E(u) 
becomes 



, ,_(n + 2M -1)! 1 

E w>(«r ^,^„ J 7 , - ^ B.i6) 



U(u) — > ,,,, , 

>E(u). 



n!(2M-l)! ((_"_),) = 



(B.17) 



This is our final result for a strict upper bound on the 
total error. Notice that this final step is equivalent to 
assuming that the contributions from all modes are as 
strong as the strongest, and add constructively. As such, 
U(u) is clearly an upperbound on the error. For one 
mode, this last step does not lead to that great of an over- 
estimate. With the addition of more modes, however, 
this step overemphasises the contribution from weaker 
modes, and leads to a (potentially much) larger overes- 
timate of the error. 
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C. THEORETICAL PREDICTION FOR N MAX 

The goal of this section is to determine the value 
of n max necessary to obtain a precision in the uv-planc 
equal to P if given u and the ionospheric modes distort- 
ing the sky. Based upon the terrible consequences which 
result from undercorrecting in the wu-plane (see section 
4), we begin with the expression just derived for the up- 
perbound on the error in hopes of avoiding this pitfall. 

In order to make sense of our expression for an up- 
per bound U(u) (equation B.16) and derive from it the 
optimal choice of n maxi we must make a few further ap- 
proximations. Some of these approximations will actu- 
ally slightly decrease the expression for the error relative 
to U(u), but are necessary in order to make sense of this 
ugly expression. 

To begin, we define 



so that 



,->\n ( n 



2M-l)\ 1 



n\{2M-l)\ ((_«_), 



. 2M 



(C.l) 



U(u)= J2 9n- (C.2) 
Next, we calculate the ratio g n +i/g n and find that 



9n+i _,(n + l + 2M-l) (2M) 
— v(u) 



9n 



(n + 1) 



(rt + 1) \n + l 



(C.3) 

where we have employed Stirling's Approximation, 



nl » (27rn) 1 / 2 n"e-™. 



Stirling's approximation is best suited for large n, but is 
actually quite accurate for small n as well. It gives an an- 
swer within 8% of the actual value for n = 1, and within 
1% for n = 9. In other words, by using this approxima- 
tion we greatly simplify our expression and sacrifice only 
a little in terms of accuracy. Using the fact that 



lim (1 + -)" = e, 

n—*oo fi 

we see that for large n this ratio reduces to 



lim 

n— »oo 



gn+i = 2Mv{u) 
g n n + 1 



(C.5) 



(C.6) 



This expression shows that at a critical value of n, 
namely, 

n c = 2Mv(u), (C.7) 

this ratio is approximately equal to 1, and g Uc +i ~ Qn a ■ 
For n < n c , g n +i > g n ] and for n > n c , g n +i < g n - In 
other words, g n is an increasing function of n until n c , 
and then decreases from then on. 

Strictly speaking, these results are only valid for large 
n >> 2M. However, in order to get an approximate 
expression for n maxi we now extend these results to all 
n. The hope is that the approximate value of n max is 
varied only slightly by this extension to small n. But 
even if this approximation causes large enough error to 
raise doubts about our quantitative results for n max , it 
should still be good enough to learn something about the 
qualitative behavior of n max . Recall that setting n max = 



is the same thing as not correcting in the uw-plane (see 
equation 2.8); 



= : 



I(u) = I(u). 



{at 



Put differently, the n = term in the sum is of the order 
of the uncorrected ww-plane, go = I{u). According to the 
above, successive corrections g n differ in magnitude from 
the previous term by a factor of 2Mv(u)/n. Therefore, 
the approximate magnitude of the term g n is 



l9n|~|/(«)| 



11/(3)1 



2Mv(u) \ (2Mv(u) 
T 

(2Mv{u)) n 



2Mv(u) 



(C.9) 
(C.10) 



We expect the distortions created by the sky to al- 
ter the magnitude of the intensity only very slightly, 
so that \I(u)\ « \Iactuai\- Furthermore, for n > n c 
(which is the case for n max ) we approximate that the 
ratio g-n+i/gn falls quick enough that we may approxi- 
mate the total remaining error as being enirely due to 
g n , g n « \I n - Iactuai\, where /„ is the intensity in the 
uw-plane after being corrected to n orders. Therefore, in 
order to obtain an fractional error / less than P for our 
uw-plane correction, we must correct to enough orders n 
so that 



\In-Iactual\ (2M V (u)Y 



tactual 



nl 



(C.11) 



(C.4) Therefore the optimal value of n max is given by 



n ma x = niin { n } such that Pn\ > (2Mv(u)) n . 

(C.12) 

Some values for n max given 2Mv(u) and P are given in 
the table embedded within the main text, table 5.1. 



D. DERIVING 10) FOR AN ALTERNATE CHOICE FOR Sd' 
A bulk of this paper has assumed that 59 takes the 



form 



M 



59' 



-ibn 



(D.l) 



This is, of course, only one of many possible choices. In 
this section, we analyze the results of instead choosing 



K 



M 



5e' = J2a k 9' x k + J2 b ™0y m > (°- 2 ) 

fe=l m=l 

where K and M are the number of terms in the 9 X and 
9 y directions, respectively, necessary to accurately model 
the distortion by the ionosphere. The a k and b m here are 
real. Similar to last time, this choice is chosen because 
it allows an analytic solution to I(u). With this choice, 
our earlier equation for I(u), 

I(u) = J ^ J d 2 e'e iff '< a '-^e- iS - 5ff 'l(u'), (D.3) 
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becomes 
I 



(5) = J J d 2 e > e ie'-(u'-u) (e-iu x Y. k * k e'*^ (D4) 
X (e-^yT, m b m e™^ (D 5) 

= j^kj ^ e <«-^) (e~i^J2 k a k e^ Iy (D6) 



But, 



Thus, 



—I{u x ,\u y ) 



= Uy 



A=l 



dut 



(D.21) 



where 



/A/' /" / \ oo / M 



(D.22) 



(D.7) 

First, we focus on evaluating I y . To do this, we first 
Taylor-expand the second exponential, 



oo / M 

n=0 \ m— 1 



(D.8) 



oo / - \ n ' M (h Q f m\l m 



,, I J. ± ; I 

n=0 l m m=l 



where the sum over l m is a restricted sum such that 
E m = Plugging this into our expression for 7 a 
gives 



x n 



'm! 



(D.10) 
(D.ll) 



: EE f ^r £/ («) /" d0; e i9 iK-^>(-™ a )X£m -^.12) 



x n (D 13) 

m=l 

(D.14) 

_ where we have assumed that the summations and the 
integrals may be freely interchanged. To solve this inte- 
gral, we first introduce an additional parameter A (which 
we will eventually set to 1) and notice that 



( D . 15 ) 



J y y J s H»,)'aA' 

= (^y^<- Xu ^ (D - 17) 
Therefore, if we define t = ^2 m ml m , then 

OO I p 

J « = ££ / du'yliu^-iuy) 71 " 1 

n i J 



n=0 lm 



X II ¥^1^ 



M 



(D.18) 
(D.19) 



=££(-*-)"-* n 



m— 1 



(b m ) lm 0* 



Plugging this into our earlier expression for we 
obtain 



M 



n— Z„ 



m=l m ' W 



where 



~27T 



(D.24) 

But / x here is of the same form as I y earlier, and there- 
fore 



oo oo / 



= EEEE-t-™*)"*"'*^-^)"*"^ - 25 ) 

n-a; n y / m Z fc 



x II -ft- II ,,, ^ z^ 7 ^ 



) (D.26) 



where 



f 'y = £ m ^ 



(D.27) 



t x — ^ ' A:/;; 



(D.28) 



and the sums over l m and h, are restricted so that 



I(u x , Aw y )(D.20) Efe h = n x and ^ m lm = %• 
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